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Introduction 


During the first year of the research program efforts have concentrated on the devel- 
opment of the numerical methods that will form the computational part of the turbulence 
closure scheme. A wave model has been developed for the two-dimensional shear layer. 
This configuration is being used as a test case for the closure schemes. Various numerical 
schemes have been examined to give efficient solutions of the Rayleigh equation for this 
geometry. These include both spectral and finite difference methods. Secondly, numerical 
methods are under development to solve the non-separable Rayleigh equation. This solu- 
tion is required for the closure scheme in more complex geometries. A model problem has 
been used to assist in the algorithm development. Two-dimensional spectral methods and 
a hybrid spectral/finite difference technique have been developed. An analytic solution 
of the Rayleigh equation for a basic elliptic flow has been obtained. This will be used to 
verify the stability codes developed for arbitrary geometries. Other numerical methods for 
solving the Rayleigh equation based on the boundary element technique have also been 
examined. These solutions are forming the basis of a model for the shock structure in jets 
of arbitrary geometry. These activities are described briefly below. 


Turbulence Closure in a Mixing Layer 

A turbulence closure scheme using a wave model is being developed. An incompressible 
two-dimensional free mixing layer has been chosen as a test case. In the future the model 
will be used to predict the mean velocity and temperature fields in jets and describe the 
characteristic frequency and wavelength properties of the fluctuating flow field. 

Since the large-scale coherent structures develop spatially, they may be best repre- 
sented by spatial instability waves characterized by the solution of the Rayleigh equation. 
At present the Reynolds stress associated with the dominant large-scale structure is mod- 
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eled by the spatially-developing most unstable mode. Future extensions of this model 
involve the inclusion of a wider frequency spectrum. 

Both the traditional shooting technique and the methods proposed by Bridges and 
Morris (1984) have been used in the solution of the Rayleigh equation. It has been found 
that a Chebyshev polynomial of relatively high order is needed to obtain accurate solutions 
if the tau method is used. The matrix operations involved in this technique, which is the 
discretization method used by Bridges and Morris, consume a relatively large amount of 
computer time. This is particularly undesirable as this solution is only one part of the 
closure scheme. Different ways to construct the coefficient matrices have been investigated. 
These include: 

(1) a finite difference method and 

(2) a Chebyshev collocation method. 

Eigenvalues of higher accuracy than those obtained using the tau method have been ob- 
tained with the order of the matrices less than half of those in the tau method. The results 
for a mixing layer with a hyperbolic tangent mean velocity profile are shown in Table 1. 
The eigenfunctions are relatively insensitive to the precision with which the eigenvalues 
are obtained so that it is unnecessary to obtain highly accurate eigenvalues. It would also 
decrease the efficiency of the turbulence closure scheme. 

The wave model has been applied locally t,o a mixing layer with the mean velocity 
profile obtained from an analytic curve fit to measured data by Patel (1973). The calculated 
Reynold stresses are shown in Fig. 1. Note that the maximum value assigned to the 
instability wave stresses is based on the measured data and neglects the contributions 
from the small-scale turbulence. The distribution shows that there are regions of negative 
turbulent production associated with the large-scale structures. This was observed by 
Komori and Ueda(1985). Moreover, it is clear from Fig. 1 that the contributions to the 
Reynolds stress from the small-scale turbulent motion must be accounted for in the model. 
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Figure 1. Calculation of Reynolds Stress Using Different Numerical Schemes 




The inclusion of the interaction between the the large-scale and the small-scale turbulence 
and the contributions of the latter to the flow development is clearly a requirement for any 
viable turbulence model. 


Method 

N 

Frequency 

Eigenvalue 

Tau 

26 

0.2 

(0.367, -0.219) 

Finite Difference 

11 

0.2 

(0.383, -0.221) 

Collocation 

11 

0.2 

(0.389, -0.221) 

Michalke(1966) 


0.2 

(0.382, -0.227) 


Table 1. Calculation of Eigenvalues 


It is recognized that turbulence comprises fluctuating motions with a wide spectrum 
of length and time scales and that different turbulent interactions are associated with 
different parts of the spectrum. Following the general ideas proposed by Hanjalic et al 
(1980) we assume that the turbulent kinetic energy spectrum may be divided into parts 
associated with the large-scale structures and that associated with the residual motions. 
Since the large-scale structures, based on the present wave model, dominate the dynamics 
mechanism of a free shear flow, we further assume that the large-scale structures are most 
effective in extracting energy from the mean flow. The turbulent energy is then transported 
to the other scales of motion through vortex stretching or higher-order instabilities and is 
then dissipated by viscous effects. A spectrum division similar to Hanjalic et al may be 
proposed as sketched in Fig. 2. ki, k 3 denote the turbulent kinetic energies associated with 
the large- and small-scale range and e; denotes the kinetic energy transfer rate across the 
inertial subrange, e represents the rate of energy dissipation in the small-scale structures. 
The transport equations for the kinetic energies may be written as 
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Figure 2. Sketch of turbulent energy spectrum 

An eddy viscosity hypothesis may be used to model the contribution to the Reynolds stress 
by the small-scale motions. 

UjUy — C\Ul 3 S{j, 

where u and l 3 are the velocity and length scales associated with the smaii-scaie motion. 
Additional assumptions are made for the small length and velocity scales, 

u - fc] /2 , 

l 3 ~ u t ~ kg^ 2 /*, 

where r is the time scale for the small-scale motions. Equation (2) can then be used to 
solve for k 3 , providing that the large-scale characteristics have been obtained from the 
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wave model and the rate of viscous dissipation is related to the rate of large- scale energy 
production by e ~ k /l. Note that the latter assumption is equivalent to 

l s ~ / ( k 3 /ki ) 3//2 . 

Extensions to this rather simple formalism are being examined. These include the 
extension of the Algebraic Reynolds Stress Closure to a multiple-scale decomposition of 
the turbulence. This extended model and the simpler model described here will next be 
applied to a calculation of the development of the two-dimensional mixing layer. 


2. Solution of the Non-Separable Rayleigh Equation 


This phase of research has centered on developing new, and improving existing, nu- 
merical methods to solve the non-separable Rayleigh equation. 

2.1 Review of Initial Progress - The first step in the solution of the Reynolds- 
averaged compressible equations of motion for jets of arbitrary geometry using a wave 
model is the description of the hydrodynamic stability of such flows. This requires the 
solution of a non-separable form of Rayleigh equation. The most unstable eigensolutions 
may then be used to model the Reynolds stress associated with the large-scale structures. 

A method has been developed to determine the eigensolutions of the Rayleigh equation 
in flows of arbitrary geometry. The equation to be solved is: 

(A - a 2 )p + (2a/n)VW ■ Vp = 0 ( 3 ) 

with boundary conditions: 

p is finite and p — ► 0 at infinity 
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where p is the pressure fluctuation, W (x, y) is the axial mean velocity, a is the axial 
wavenumber (the complex eigenvalue), u> is the wave frequency, and Q = w - oW. 

In order to test various numerical algorithms for solving this problem a model problem 
with a known analytic solution has been posed. The boundary value problem is given by: 

A <f> — 2 au>(d x (f) + d y (f>) — 2a 2 <t> = 0 (4) 


with, 


<f> = 0 on dtt 


where 

H:= {(x,y) € R 2 |-1 < x,y < -1} 

The method that has been used to determine the eigenvalues a is a spectral method using 
a two-dimensional Chebyshev series approximation of the form: 


EE < fimnTn ( 5 ) 

Here the summations are taken from 0 to N. 


The determination of a goes as follows: i) The P.D.E. is integrated to eliminate all 
derivatives of the dependent variable, ii) This integral equation is discretized using the 
Chebyshev series, equation (5). iii) The boundary conditions are discretized, iv) Steps ii) 
and iii) reduce to the problem of determining all a such that: 

det[Aa 2 + Bo + C] = 0 (6) 

where A, B, and C are (N + l) 2 x (N + l) 2 matrices depending on the discretization, v) 
Equation (6) is solved using globally or locally convergent schemes. 

2.2 Review of Current Research - Logically, the next step is to generalize the 
Chebyshev integral method from the model problem (4) to the Rayleigh equation (3). 
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However, this transition is complicated by the non-constant coefficients present in (3). To 
realize the problem, notice that in general the variable coefficients of the Rayleigh equation 
must be expanded in infinite series of the form (5). Then the coefficient series must 
be multiplied by the respective series approximating a term in the Rayleigh equation. 
Algebraically, the multiplication of multiple infinite series is very complicated and the 
resulting formulas are tedious to code. Therefore, a spectral method is sought to solve the 
Rayleigh equation which avoids the algebraic complications of multiplying multiple infinite 
series together. Such a method is outlined below. 

Another difficulty associated with the solution scheme developed for the model prob- 
lem (4) is the size of the resulting lamba matrix, (6). It is anticipated that a typical value 
of N will lead to a very large matrix. Hence, from a computational viewpoint it would 
be advantageous to reduce the order of the lamba matrix. Currently, methods combining 
finite differences and pseudospectral approximations are being developed. These methods 
have the potential to reduce the order of the lamba matrix from ( TV -I- l) 2 x (TV + l) 2 to 
(TV + 1) x (TV + 1). Schemes combining finite differences and spectral approximations, i.e., 
hybrid schemes, will be discussed in section 2.2.2. 

2.2.1 Pseudospectral Chebyshev Methods - Let L denote the partial differential oper- 
ator associated with the model problem (4). A scheme is sought to determine the complex 
wave numbers a of the model problem such that 

L<f> N (xij,u,a) = 0, (7) 

is satisfied exactly, where <T>jv is a pseudospectral approximation of the solution. Here, the 
frequency w is fixed, and x;y are collocation points in fi for 0 < i,j < TV defined by: 

X t y = (x,,yy), 

= (cos[7ri/TV],cos[7ry/TV]). 
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Furthermore, <f>u is defined by, 


<f>N{x,y) = £X>(x«)/, («)/,(*). ( 8 ) 

where fj(z) (z = x or y) is given by, 

fM = ((i - ^)r;w(-i)' +I )/(c^ 2 (z - *>)), (9) 

with c 0 = c N = 2, and cy = 1 for j = 1, • • • , iV — 1. 

The discretization of eqn. (4) is then given by evaluating eqn. (7) at the collocation 
points in the interior of Q and combining these equations with eqn. (8) satisfied at the 
collocation points on the boundary of fi. The advantage of the pseudospectral Chebyshev 
method becomes clear when the operator L is defined to be the Rayleigh operator. Instead 
of expanding the variable coefficients of L as an infinite series, the coefficients are simply 
evaluated at the collocation points. Thus, the algebraic problem encountered in the integral 
method is avoided. 

2.2.2 Hybrid Schemes - In order to reduce the order of the discretization matrix, 
resulting from a multiple series expansion, one variable expansions are currently being 
used to develop a method to solve the model problem. In general, a one variable series 
will generate a smaller discretization matrix. 

Consider a solution of the form: 


<l>N{x,y ) = ]Pa ( x t ,y)/ t (x). (10) 

Here the summation is taken from 0 to N. Moreover, the /,•’ s are functions which may be 
used to approximate <f>. 

To outline a hybrid scheme, let the functions, /,-, in eqn. (10) be defined by eqn. 
(9). Therefore, eqn. (10) provides a pseudospectral approximation to the solution in the 
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x direction. In the y direction, a finite difference scheme will be applied. Let V 1 (y) and 
V 2 (y) denote arbitrary finite difference operators in y for the first and second derivatives 
respectively. Then substituting (10) into the model problem and applying the V l (y)’s 
gives, 

[V 2 (y) - 2waV 1 (y)] <f> N = [2wad« + 2 a 2 - d xx ] <j> N , (11) 

where the operator on the right hand side of eqn. (11) is approximated by the psuedospec- 
tral method. 

The boundary conditions are established by recasting the boundary value problem 
in y as an initial value problem. Since the derivative of <t> with respect to y is unknown 
along y = — 1 the solution is constructed as a sum of linear problems for which the initial 
conditions are 

The solution on y = 1 at the collocation point X{ is given by 

<t> N {xi,l) = '^T,A j <i> , N {xi,l). (12) 

Application of the boundary condition <j>(x, 1) = 0 yields a system of homogeneous simul- 
taneous equations for the Aj. For a non-trivial solution the determinant of the coefficient 
matrix given by these equations must be zero. Only when a is an eigenvalue will this 
condition be satisfied. 

This method is readily extended to the boundary conditions posed for the Rayleigh 
equation. Casting the non-separable boundary value problem in this way should reduce 
the operation count by a factor of N 2 and yield a much more efficient evaluation of the 
eigenvalues. 


3. Boundary Element Methods and Shock Structure Tam and his co-workers have 
shown that both the large-scale coherent structures and the shock structure of jets may be 
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modeled using a wave analysis. The essential difference between the two cases is that the 
former are modeled as travelling waves while the latter are standing waves. In this section 
a numerical method, the boundary element method, has been used to determine the shock 
structure in jets of arbitrary cross-section modeled by a vortex sheet. In the next stage 
of the analysis the effect of finite mixing layer thickness and its axial variation will be 
considered. Modeling this physical problem provides a vehicle to test different numerical 
methods for the solution of the Rayleigh equation than those described in Section 2 as well 
as giving additional insight into the physical structure of jets of arbitrary shape. 


The shock-structure is to be analyzed for multiple jets of arbitrary shape. The model 
developed by Tam and his co-workers, e.g. Tam et al (1985), is used here. This model 
incorporates the effects of a finite thickness mixing layer and the axial flow divergence. 
We are concerned with solving the equations of linear hydrodynamic stability in which the 
mean flow characterstics are arbitrary functions in a plane normal to the jet axis or axes 
and a slowly-varying function of the axial distance. A typical cross-section of a single jet is 
shown in Fig.3. This defines three regions. In regions I and III the mean flow properties are 


constant. These regions correspond to the potential core of the jet and the ambient fluid 
surrounding the jet. Region II represents the jet shear layer in which the mean velocity 



Region II 

(variable 

coefficients) 


Region III 

(constant 

coefficients) 


Fig. 3: Sketch of Regions of Jet Cross-Section 
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The equation for the pressure fluctuations to be solved in these regions is given by, 


Ap + 


2k 

{P'kU) 


Wdp dU_dp 
dy dy dx dx 


- [k 2 — M 2 (f3 — kU) 2 } p = 0, 


(13) 


where, p is the pressure fluctuation, k is the axial wavenumber, (3 is the frequency, and 
U(x,y) is the mean axial velocity. Equation (13) is obtained by taking the divergence of 
the linearized momentum equation and using the continuity and energy equations. 


The analysis of the problem is carried out in the following stages: 

(i) A linear shock cell model is used as a first order approximation to analyze regions 
I and III. In this model, the mixing layer of the jet, which is assumed to be thin, 
is approximated by a vortex sheet. This assumption eliminates the need to analyze 
region II. The details of this model and the solution technique are discussed below. 

(ii) The inclusion of the effects of the finite thickness shear layer and that of axial flow 
divergence. The inclusion of these effects require a numerical solution in region II. 
The solution obtained for region II has to be matched with the solutions of regions I 
and III at the common boundaries. 

(iii) Extension of the model to include the effects of multiple jets of arbitrary cross-sections. 


3.1 Vortex Sheet Shock Cell Model - Consider a shock cell system in a jet column 
bounded by a vortex sheet as shown in Fig. 4. For convenience, a Cartesian coordinate 
system centered at the nozzle exit with the x-axis in the direction of the jet centerline 
will be used. The surface of the vortex sheet bounding the fully-expanded jet is given by 
<S'o(y,«) = 0. There is no disturbance outside the jet. The linearized equations of motion 
inside the vortex sheet are: 


dp 

„.v. v + «,-=°, 

dV 

= - Vp ’ 


P = a 2 p. 


(14) 

(15) 

(16) 
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pj, Uj and ay are the density, velocity and the speed of sound of the fully expanded jet. p, p 
and V are the density, pressure and velocity associated with the linear shock cell structure. 



Fig. 4: Vortex Sheet Shock Cell Model 


From eqns. (14) to (16) it is found that the pressure p satisfies the equation 

Vp 




( 17 ) 


with p = 0 on the boundary So(y,z) = 0 and at x = 0, p = A p inside So{y,z) = 0 and 
V ± = 0. 


A general solution of the vortex sheet shock cell boundary value problem can be found 
by writing the pressure fluctuation as: 

p{x, y, z) = 4>{y, z) cos(fcx). (18) 

where k is an as yet unknown axial wavenumber. The equation for <j>(y, z) may then be 
written 

V<j>+\ 2 <f) = 0 where, \ 2 = (M 2 — l)k 2 , (19) 

with <t> = 0 on So(y,z) = 0. This is an eigenvalue problem with k as the eigenvalue, which 
is to be determined. 
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The Boundary Element Method is used to solve the Helmholtz equation (19). 
This integral technique is particularly suited to problems posed in an arbitrary geometry. 
The development of the method is discussed below. 

Consider an arbitrary domain as shown in Fig. 5 where the boundary is divided into 
N panels. £ = £(ai,<Z2) are the locations of the node points and a = 0(61,62) are the 
locations of the mid points of each panel. 




Fig. 5: Sketch of Boundary Divided into Elements 

Let F(x | y) be the fundamental solution of the Helmholtz equation, i.e. 

(A + A 2 )F(x|y) = -6(x-y), (20) 

where, x = (xi,X2) and y = (s/i , j/2) and x and y are arbitrary. Then, 

F(x\y)= jff ( < 1 , (A|*-s,|) 

Application of the divergence theorem and noting that <f> satisfies eqn. (16) and is zero on 
the boundary yields, 

[ ~r^F(a \£)do for * = 1,2,3, ,jV 

2 Js„= 0 ov 
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Let Li be the length of the panel on the arc T t . Then, 
for t=j, 



For i^j. 


<“>(*) = j {f K’W* - 01 ) +4ff^ ) (A|a, - a, |) + A|a, - £ y+1 |)] j , 
where, H(z) are Struve functions that can be expressed as a series of Bessel functions. 


The eigenvalues are obtained by a local iterative scehme which is given by, 

A*+i = A k - l//(Afc), 


where, 

/(A fc ) = Tr[ M - 1 (A fc )/(A,)], 

Tr denotes the trace of a matrix and a prime denotes the derivative of the matrix with 
respect to A*. 
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Thus the eigenvalues can be determined starting from an initial guess. 


The Boundary Element Method has been used to solve the vortex sheet shock cell 
model for rectangular and circular jets. These two geometries are chosen for the reason 
that their analytical solutions for the eigenvalues are easily obtainable. The numerical and 
the analytical values are compared for the smallest eigenvalue as this primarily determines 
the shock cell spacing. The results are shown in Tables II and III. 

As is to be expected the number of panels has a great effect on the accuracy and the 
number of iterations required. A compromise between the accuracy required and the total 
CPU time has to be reached. 


Analytical value of A = 4.442883. 


Number of 

Number of 

^ numerical 

^ initial 

% error 

Panels 

Iterations 




4 

29 

(5.14421,0.1056) 

(4.0, 0.0) 

15.78 

8 

11 

(4.31998,-0.05169) 

(4.0, 0.0) 

2.77 

12 

11 

(4.40674,-0.02355) 

(5.0, 0.0) 

0.81 

20 

11 

(4.43237,-0.00627) 

(5.0, 0.0) 

0.24 

28 

12 

(4.43795,-0.00243) 

(5.0, 0.0) 

0.11 

40 

13 

(4.44057,-0.00086) 

(5.0, 0.0) 

0.05 


Table II: Calculations for a Rectangular Jet 
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Analytical value of A = 2.40482. 


Number of 

Number of 

^numerical 

^initial 

% error 

Panels 

Iterations 




4 

28 

(3.63751,0.07463) 

(3.0, 0.0) 

51.25 

8 

12 

(2.58661,-0.02283) 

(3.0, 0.0) 

7.60 

12 

12 

(2.47651,-0.00827) 

(3.25,0.0) 

3.00 

20 

11 

(2.42874,-0.00201) 

(2.75,0.0) 

1.00 

30 

12 

(2.41510,-0.00062) 

(2.75,0.0) 

0.40 

40 

14 

(2.41000,-0.00027) 

(2.75,0.0) 

0.20 


Table III: Calculations for a Circular Jet 


In the next stage of the analysis the effect of the finite thickness of the mixing region 
will be examined. 
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